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Abstract 


For 3D vector fields, the governing formula of invariant manifolds grown from a hyperbolic cycle 
is given in cylindrical coordinates. The initial growth directions depend on the Jacobian matrices of 
Poincaré map on that cycle, for which an evolution formula is deduced to reveal the relationship among 
Jacobians of different Poincaré sections. The evolution formula also applies to cycles in arbitrary finite 
n-dim autonomous continuous-time dynamical systems. Non-Möbiusian/Möbiusian saddle cycles and a 
dummy X-cycle are constructed analytically as demonstration. A real-world numeric example of analyz- 
ing a magnetic field timeslice on EAST is presented. 
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1. INTRODUCTION 


In the tokamak community, the magnetic topology is, most of the time, assumed to be nested flux 
surfaces, e.g. in Grad-Shafranov equation, EFIT, VMEC, etc. Based on this assumption, people 
have established a dedicated theory [1] of magnetic coordinates in which all magnetic field lines 
on a flux surface are straight. Nevertheless, the three-dimensional (3D) effect is ubiquitous in 
real-world fusion experiments since any facility in fusion machines, except the central solenoid 
and poloidal field coils, is 3D. For example, the toroidal field coils exhibit a ripple effect, the mi- 
crowave and radio-frequency wave heating impose their distinctive localized distribution pattern 
of the induced current, and the neutral beam injection (NBI) has an obvious non-axisymmetric 
current effect. Furthermore, most instability modes of plasma, such as the tearing mode and 
sawtooth mode, imply a significant 3D topology change of the magnetic field. The 3D effect is 
unavoidable in the magnetically confined fusion research, thereby requiring a deeper comprehen- 
sion on the global field structure. 

An enormous amount of fusion research has attempted to stimulate a chaotic field layer at 
the plasma boundary by either resonant magnetic perturbation (RMP) coils [2-4] or other means 
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[5] to mitigate destructive type-I edge localized modes (ELMs). The theoretical basis was estab- 
lished in 1983 by Cary and Littlejohn [6] to estimate how wide the island chains are when an 
axisymmetric magnetic field is perturbed by a non-axisymmetric one, after which hundreds of 
researchers implement RMP coils to mitigate and suppress ELMs [7,8]. The method is called mag- 
netic spectrum analysis nowadays, heavily relying on Fourier transform of the radial component of 
the perturbation field, which is a linear operation in functional spaces. Therefore the utility of 
magnetic spectrum analysis is limited inside the plasma and becomes less and less accurate as 
the perturbation is strengthened, because Fourier transform is merely a linear operation and not 
capable of explaining the nonlinear behaviour. 

With the aid of modern dynamical system theory, the structure of a 3D vector field can be 
comprehended and analyzed in terms of invariant manifolds [9,10]. Various numerical methods 
have been developed to grow them [11-16]. Kuznetsov and Meijer systematically investigated 
the bifurcation behaviour of 1D and 2D maps when their parameters change [17, 18], presenting 
diverse analytic and numerical methods to study maps. However, the methods taken by math- 
ematicians are too general to capture the essence of a 3D vector field, leaving the analysis as 
complicated as before. 

Since the magnetic field dominates the plasma transport in magnetically confined fusion ma- 
chines, the evidence of transversely intersecting invariant manifolds shall be easy to observe, 
which is a signature of chaos. The invariant manifolds are essential to determine the chaotic field 
regions, which induce a mixing effect inside the plasma. The plasma edge is not suitable to be 
characterized by a single closed surface when the 3D effect is strong, for which it is proposed 
to use the notion of invariant manifolds of outmost saddle cycle(s). In fact, it has been observed in 
some existing simulation [19-22] and experiment results. On EAST, the helical current filaments 
induced by the lower hybrid wave heating impacted the plasma edge topology and caused an 
evident splitting of strike points in experiments [5]. 

If RMP coils are imposed to suppress ELMs in tokamaks, the heat flux pattern at the diver- 
tor exhibits a complex toroidally asymmetric distribution, which poses challenges to ITER and 
DEMO divertor designs [23]. Simulation results demonstrate that the divertor plasma regions 
with connection to the bulk plasma are dragged further outside when the asymmetry is inten- 
sified [24,25]. The field line connection length and the magnetic footprint (how deep the field 
lines penetrate into the bulk plasma) distribution at the divertor are usually ribbon-like. The 
RMP does mitigate or even suppress the ELMs, otherwise which can cause intolerable transient 
particle and heat flux. On the other hand, significant heat fluxes may arise far from the strike 
point originally designed for axisymmetric cases [26], possibly damaging the fragile parts of the 
divertor. Moreover, it remains unknown whether the total heat flux leaked from the bulk plasma 
is reinforced by the perturbation. 

Previous research has attempted to draw the two transversely intersecting manifolds of hy- 
perbolic cycles for 3D toroidal vector fields. Ottino from the community of fluid mechanics [27], 
Roeder, Rapoport, and Evans from the fusion community [28,29] have drawn up the relevant 
figures. Abdullaev has deduced an approximate (to first order in € because Poincaré integral 
is used) analytic implicit expression of the invariant manifolds of the outmost X-cycle when a 
single-null configuration is perturbed [30]. This paper carries forward the research and directly 
deduce the intrinsic analytic formula of the invariant manifolds of hyperbolic cycles without need 
for approximation. 

Sect. 2 explains the definitions and denotations used in this paper. The theory summary 
of this paper is put in Sect. 3, while the detailed proofs are put in Appendix A. Sect. 4 offers 
intuitive examples to help understand the theory. Non-Möbiusian/ Möbiusian saddle cycles and 
a dummy X-cycle model are constructed analytically. A real-world numeric example of analyzing 
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a magnetic field timeslice on EAST is presented subsequently. Sect. 5.1 and Appendix B give some 
comparisons with others' works. Sect. 5.2 is the conclusion. 


2. DEFINITIONS AND DENOTATIONS EXPLAINED 


Vector fields are presumed to be at least once continuously differentiable, i.e. of class Cl. Al- 
though a vector field is denoted by B, it is worth emphasizing that it does not need to be 
divergence-free in this paper. 

The flow X(xo, f) induced by the field B and the corresponding flow X jj (xo, Ps, fe) in cylin- 
drical coordinates play a great role in the deduction of theory, especially the latter one. (A flow 
is often denoted by déi or $(t, xo) in other literature, but p has been used as the azimuthal angle 
in this paper. X and X,,; denote the field line tracing (FLT) flows instead.) Denotations used are 
listed in Table 1. The symbol D serves as a differentiation operator w.r.t. the initial condition x9, 
as distinguished from those w.r.t. time-like variables t, ds or Qe. 


Table 1: Denotations 


Cartesian Cylindrical 
E dxr/dx =  RBg/B 
FLTODE system a pa [17 Rus i i () 
z= Z dxz/dxg = RBz/Bg 
> 
Init condition xo = (Xox, Xoy, Xoz) xo = (Xon, Xoz) 
Flow X(xo, t) X pol (X0, $s, Pe), subscripts s for start, e for end 
RB 
dX (xo, t) = B(X (xo, t) da, X pol (x0, ds, Pe) = x (X (xo, s, Pe) 
$ 
Differentiation 9X(xg, f) IX poi (Xo, Ps, Pe) 
DX(xo, t) = —————. DX REL 
wrt. xo (xo, £) O(Xox, Xoy, Xoz) poo, Ps, Pe) O(XoR, Xoz) 
d d 
By chain rule 3; PX Gro, f) = VB(X(xo, OU DX(xo, t). ae P X vot (0 Ps, Pe) = A(ĝe) DX poi (x0, Ps, Pe), 
(2) (3) 
A (RB,u/B 
pol/ 24 
where A(fe) := ) (X poi (Xo, Ps, Pe), Pe) 


a(R, Z) 


Naturally, the Poincaré map P(x0, $), xo = (Xor, xoz), is defined by the standard R-Z semi- 
infinite planes Zi) at $ angles, recording the strike points crossing the planes in one direction. If 
a trajectory flies around Zi) through the opposite semi-infinite plane X(¢ + zt), then by definition 
the Poincaré map does not record the strike point on Zb + 7). Since the flow X ,jj(xo, ps, Pe) and 
the map 7 (xo, di are used frequently in this paper, different terminologies are used to distinguish 
the continuous-time dynamics from the discrete-time one, e.g. equilibrium point, trajectory and cycle 
are for continuous-time dynamics, while fixed point, orbit and periodic orbit are for discrete-time 
one. 

This paper follows the definition of stable and unstable manifolds given by Jacob Palis and 
Welington de Melo [31], p. 73 for a hyperbolic fixed point (p. 59) and p. 98 for a hyperbolic 
cycle (p. 95). In both cases the manifolds are defined by w- and a-limit sets, so one can unify 
the manifold definitions in the continuous and discrete dynamics by considering manifolds of 
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a hyperbolic invariant set S (a def. at p.160 for the case of a map), since both hyperbolic fixed 
points and hyperbolic cycles are hyperbolic invariant sets. Note a hyperbolic periodic orbit is 
also a hyperbolic invariant set, which can be handled by the unified definition. Let (M, f) be 
a dynamical system, where f is a map P or a vector field B defined on the manifold M. The 
stable and unstable manifolds of a hyperbolic invariant set 5 for f are defined by w- and a-limit sets 
respectively 


W*(f,S):- (p € M : olf, p) = S, (4) 
WES) = ip € M : a(f,p) = Sj. (5) 


The term hyperbolic ensures the nearby trajectories/orbits on stable (resp. unstable) manifolds 
approach (resp. get away from) 5 at an exponential rate, without which the sets defined above 
are only qualified to be named stable (resp. unstable) sets [32] p. 257, because the stable manifold 
theorem ([31] p. 75) requires 5 to be hyperbolic. Colloquially, the stable (resp. unstable) manifold 
is the set of points that will flow into (resp. flowed out of) S at an exponential rate; see [32] p. 258. 
Kuznetsov and Meijer simply write JE) > S as k — oo (see [18] p. 4, f is a map) without 
explaining what they mean by "converging to a set S", which is suspected to be an informal 
denotation of w- and a-limit sets. 

For notational convenience, arguments can be omitted if suitable, e.g. DX(f) and DX poi($s, Pe) 
are short for DX (xo, dr) and DX poi(xo, Ps, Pe) by omitting xo, DP" (p) for D'P" (xo, $), and we/s(s) 
for W"/5(f, S). 


3. FROM FIELD LINE TRACING TO INVARIANT MANIFOLDS 


DX and DX ai imply the change of differential volume and area during FLT, respectively. Sup- 
pose a 2D map is written as (x, y) +> (u,v). The differential area expands, shrinks, or keeps con- 
stant after being mapped, as revealed by the following exterior product of differential 1-forms, 


du ^ dv = (ar 24 3,49) ^ (Zax i ) dx ^ dy. (6) 


As X pol(%s, Pe) is a typical 2D map from the section X(Ps) to (pe), the determinant of DX por (Ps, Pe), 
denoted by [DX poi( Ps, de) 


, is indeed the same thing as du dun — dyu 0,0. One could be curious 
about the geometric meaning of GEN and conjecture that it must be related to the 
divergence of the field, since it has been well-known (see [33] p. 408) that for |DX(xo, t)| 


|DX(xo,t)| = eho t VB(X(x0,1))dt — el VBX (x07) (7) 


which indicates that, for a divergence-free field, |DX(xo, t)| is always zero. A similar formula for 
DX yoi($s, 90) is deduced (proof in Appendix A.1) to reveal the relationship between [Dx pol (Ps, Pe) 
and the divergence along the corresponding trajectory X pol (Ps, PD), Ps € $ < Pe, as shown below, 


E be R(V - B) , \ Bolo, 
= exp ( f, ay) M (8) 


[DX por ps, de) 


which applies to 3D vector fields of class C 1, no matter whether they are divergence-free or not. 
More importantly, DP*""(xo,@) with xo on an X-cycle y of m toroidal turn(s) decides the two 
X-leg directions of the X-point. DP" (xo, di = DX pol(Xo, $, $ + 2mn) if By is positive everywhere. 


4 
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It is the two eigenvectors of DP" (xo, di) that dictate towards which direction the two invariant 
manifolds of that hyperbolic cycle grow at the beginning. However, DP=" (xo, 4) itself is more 
difficult to solve for than its determinant, which is a constant independent of ¢ for the cycle. This 
is because the right-hand side of Eq. (8) becomes a constant as fe = Ps + 2mm, i.e. 


$.--2m7t F ` 2 
DX poi($s, Ps *2mn) — exp (/, "a Pap) är. exp (f E) . (9) 


S 


Unlike the determinant |DP*"(x0, $)|, the matrix D'P^" (xo, $) varies w.r.t the azimuthal angle 
$. The evolution rule of DP^" (xo, $) w.r.t. p is revealed in section 3.1. 

To calculate DX poi(Ps, Pe) by integrating Eq. (3) requires that By on the trajectory does not 
change its sign. Otherwise, RB,,1// By would be undefined due to zero By. These trajectories are 
not useless and could be well-defined in Cartesian coordinates. Suppose a trajectory goes from 
Ps to pe, during which By may change its sign for several times. The zero By singularities cause 
inconvenience to solving for the DX ail, Pe). To get around the zero By singularity issue, one 
can solve for the corresponding DX first and then change its coordinates back to the cylindrical 
system. The following DX to DX ,,; formula (proof in Appendix A.2) tells how to do so, 


sl 


UB E $ —Rsing cos $ —Rsin j A 
sind R cos $ DX |sing Rcos $ = [DX po N | ; 
1 E RBz 1 1 
Bo end end start 
(10) 


where the subscripts start and end mean that the corresponding matrices are evaluated at the start- 
ing and ending point of this trajectory, respectively, i.e. (x pol (Ps, $s), Ps) and (x pol (Ps, Pe), fe): 


3.1 The evolution of DP~"” along a cycle 


For a cycle of m toroidal turn(s), the relationship among the DP~"(#) matrices at neighboring 
sections is desired, without which the calculation of their eigenvectors would spend unnecessarily 
enormous computational resources. The most primitive approach is definitely repeating the 
integration of Eq. (3) from Ọs to Ps + 2mn once and once again for various $s, ie. from dn 
to du + 2mm, from $» to $» +2mn... To avoid this horribly primitive approach, the following 
DP*" evolution formula is deduced (proof in Appendix A.3) to reveal how DP*""(¢) varies along 
the cycle, 


d 
ag 
where the square bracket denotes the commutator, i.e. [A, B] = AB — BA. In the dummy X-cycle 
demonstration Fig. 2(a, b), Sect. 4.2, arrows are drawn to indicate the directions of eigenvectors. 
The DP*" evolution formula can be applied to cycles in autonomous n-dim flows (n > 2 and 
finite). For other dimensions than n = 3, the denotation |DX(xo, T)| is preferred than |DP"(4)], 
where T is the period of the cycle, because the former one does not rely on the choice of Poincaré 
section. The n-dim version of Eq. (11) in Cartesian coordinates is shown below 


DP-"(9) = [A(9), DP", (11) 


E DX (X, t), T) = [V B(X(xo, t), DX(X(xo,t), T)] (12) 


Furthermore, it is desirable to deduce how an eigenvector of DP" (0) evolves along an X- 
cycle, so that one can get rid of the arbitrariness of the computed eigenvector direction, which 
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depends on the specific eigen-decomposition numeric algorithm. If the eigenvector rotates a lot 
during evolution, the employed numerical method might give a reversed direction without con- 
sistency, i.e. the computed eigenvector may jump to the opposite side suddenly. The following 
DP" eigenvector evolution formula (proof in Appendix A.4) extracts the underlying rule govern- 
ing the rotation of D'P"() eigenvectors along the cycle. Let the eigenvectors of D^" ($) be 
denoted by v; — [cos 0;($), sin di) d i € {1,2}. The derivative of O($) :— diag(01, 02) w.r.t. $ 
satisfies 


= 


e = (|, Tva- pp | "iv (Dp*") v, (13) 


1 
where V :— [v1, v2] and A := diag(A1,A2). We discovered in numeric implementation that the 
formula above (13), though accurate, but encounters some numeric issue while handling the X- 
cycles in island chains, because the two eigenvectors are so close to each other that some matrices 
in this formula might become pretty singular, i.e. have big conditional number. A much more 
robust way is to directly use DP" (#) evolution formula (11) to evolve DP”"(4) and later comb the 
directions of eigenvectors. 

Traditionally, fixed points of 2D maps are classified into hyperbolic, elliptic, and parabolic 
types based on their Jacobian eigenvalues. For the cycles of 3D flows, the authors want to imitate 
this naming convention. It is worth emphasizing that the eigenvalues of D'P^"' (i) keep constant 
during evolution.! Hence, it is safe to classify a cycle y of m toroidal turn(s) by its DP” eigen- 
values. The A-invariance ensures the safety of such classification, because one does not need to 
worry about that the DP" (o) eigenvalues at different $ are different. 

If both eigenvalues of DP” are not on the unit circle S of C, the cycle is said to be hyperbolic.” If 
only one eigenvalue on the unit circle, the cycle is called partially hyperbolic (but not hyperbolic). 
If both eigenvalues are on the unit cycle of C but neither equal 1 nor —1, the cycle is defined to 
be elliptic. If the two eigenvalues are identical to 1 or —1, the cycle is defined to be parabolic. 

Furthermore, a saddle cycle is defined to be one with |A;| < 1, |A2| > 1. Those saddle cycles 
with both eigenvalues negative (resp. positive) are called Möbiusian (resp. non-Möbiusian). Note 
that the Móbiusian cycle defined here is different from the classical Móbiusian strip. The cycles 
with both A inside (resp. outside) the unit circle S of C are defined to be sinking (resp. sourcing) 
cycles. 


1This is a well-known fact in the ODE community, because it is easy to verify that DP(p2) and DP ($i) are similar by 
checking DP(42) = DX pol ($2, 61)  DP($1)DX pol (2, $1). The commutator form of the right-hand side of the evolution 


EM 


formulas (11, 12) also ensures the A-invariance. Let x; be a right eigenvector of DP" and y the corresponding left 
eigenvector, 

DPH" Xi = AiXi, y DPH" = Au. (14) 
Given a univariate matrix function A(¢), the derivative of its eigenvalue w.r.t. the single parameter d satisfies AL = y A! x; 
[34]. Now substitute A for DP^" as shown below. It immediately lets us know both At) of D^" (i) do not change 
with $. 


A = yl (DPH"Y xi 
= yL(A DP" — DP" A)x; = yT A(A;xj) — (Aiy?) Ax; = 0. (15) 


?This definition is equivalent with the conventional one [31] p. 95, whose Poincaré section is chosen to be transversal 
to the local vector field. 
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non-Móbiusian if both A positive 
saddle ME . . 
. Móbiusian if both A negative 
hyperbolic EN 
sinking ^ if both A inside S 


cycles of sourcing if both A outside 5 
3D flows | partially hyperbolic 


guef if one eigenvalue A = 1 or —1, while the other A € R \ {1,—1} 


llipti if both A = 
ban hyperbole f iptic if bo on Sbut £ +1 


parabolic  ifbothA —1or —1 


Magnetic fields are typical divergence-free fields, in which the so-called X-cycles, O-cycles, 
and the cycles on rational flux surfaces can now be formally defined as hyperbolic (meanwhile 
saddle), elliptic, and parabolic, respectively. 


3.2 Invariant manifold growth formula in cylindrical coordinates 


Consider that an invariant manifold of a hyperbolic cycle y may grow endlessly, then one of the 
two parameters of the manifold is naturally chosen to be the arc length s of the curves intersected 
by the 2D manifold W"/*(B, y) and R-Z cross-sections. For the other coordinate, the azimuthal 
angle $ is chosen. It is defined that s — 0 on the cycle and s increases towards the positive 
infinity as the manifold grows away from the cycle. The diagram (Fig. 7) in Appendix A.5 could 
be helpful for readers to understand the geometry, which illustrates the relationship among the 
differentials used. The diagram is put in Appendix because those readers who follow the proof 
there would need it more. 

To be accurate, this paper defines a stable/unstable (manifold) branch to be a connected com- 
ponent of the manifold minus its invariant set, i.e. a connected component of W'/5(f, S) \ S. For 
example, a saddle fixed point of a 2D map has 2 eigenvectors and 4 invariant branches. An invari- 
ant branch of y is parameterized to be xus (s,$) — (xil *(s, p), x SC? 9) , Where the superscripts 
u and s indicate whether the branch is unstable or stable. To deduce the governing equation of 
X"/5(s,p), one simply needs to analyse the differential relationship appearing in FLT, which is 
concluded in the following invariant manifold growth formula (proof in Appendix A.5, requiring 
no more than knowledge of multivariable calculus), 


RB, axu/s 
" — RE (X5, — 66,9) 
e e (16) 
de , A RB, "m E axu/s (5 
= By Hu ap ‚$ A 


with the initial condition 0;X"/*(s, oi i set to be the normalized eigenvector of DP" (p). Nat- 


urally, dX "its, di is 2m 7zt-periodic in $ for a non-Móbiusian saddle cycle. The denominator is 
essentially ds/d, so the sign + shall take + if the field line moves away from the cycle as $ 
increases. Otherwise, it takes —. 

For Móbiusian saddle cycles, the invariant manifolds can be grown similarly with some subtle 
difference. A non-Möbiusian saddle cycle has two invariant branches for each (un)stable man- 
ifold. With regard to a Móbiusian saddle cycle, the two branches of a (un)stable manifold are 
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considered as a whole (since they are connected). Double the period of X "its, p) in $ from 2m7 
to Amn. X"/ *(s, 6) is opposite xul *(s, $ + 2m 7t) across the cycle y. Then the growth formula 
works again. 

The growth formula would not grow a rational ? flux surface from a parabolic cycle on that 
surface, since any field line does not cover that surface. To grow this surface, the DP*" evolution 
formula is sufficient. One simply needs to move the cycle in the directions of the DP" eigenvec- 
tors step by step. The role of D" evolution formula is to accelerate the computation of DP” (o) 
at all sections. For an irrational flux surface, one can choose a non-invariant “cycle” which does 
not obey the FLT ODE system (1) and then employ the growth formula. For example, in an ax- 
isymmetric field, pick a "cycle" whose R, Z coordinates are constants. The growth formula then 
degenerates into 


0 
oxu/s RB 


pol ox” o 
= ted, (16 revisited) 
E n l-l 


therefore 9 X''/5 /ds is parallel to B,,; and of unit length. 

If what is known is limited to one section Zi) (so $ is omitted in this paragraph), then the 2D 
invariant manifold growth formula degrades to a delay ODE describing 1D invariant manifolds 
for a 2D map. Let (xj)? , be a hyperbolic periodic orbit under P. Parameterize an 1D invariant 
branch of W"/ °(P™, x1) by its arc length s as Kei *(s) : R39 > RÊ, whose inverse is denoted by 
s(X). Then one can acquire the following equations to grow the manifold by simple analysis 
along the branch: 


dX"(s) M dx" i 
a = DP"(P-"(x'G))- 3 (s (P^" (x"(9))) / AT Ka 
dX*(s) P dX? 
d; = DÉI "UP (x*()) a (P" (x*()))) A l+ lla» SS 
where the ellipsis dots --- in the denominators, which serve to normalize, denote the same thing 


as the numerators. Both equations above also hold for an invariant circle C. The only thing that 
needs to be modified is that the circle should be parameterized as a function X(s) : IR — C C IR? 
of period the circumference of the circle, whose inverse s(X) : C — R is now a multivalued 
function. 


4. DEMONSTRATION OF CYCLES AND INVARIANT MANIFOLDS 


Having developed the systematic theory to characterize the invariant manifolds in 3D autonomous 
flows, the formulas have been implemented to present readers with a vivid picture. Two analytic 
and one real-world examples are exhibited in this section. As the simplest model, the first ana- 
lytic model is a saddle cycle as shown in Fig. 1, which can be either non-Móbiusian or Móbiusian. 


3In KAM theory, if there exists a diffeomorphism q : 7^ — T’ ( T is the standard d-torus) such that the resulting 
motion on T^ is uniform linear but not static, i.e. dg(x)/dt = w € R? \ {0} is constant, then w € IR? is called a frequency 
vector (a.k.a. rotation vector) of 7^. A frequency vector w is said to be rationally dependent or commensurable if there 
exists a vector k € Z’ \ {0} such that k -w = 0. Specifically when d = 2, this paper defines T? to be a rational (resp. 
irrational) flux surface if the corresponding w is commensurable (resp. incommensurable). The frequency vector w is 
unique up to a transform of multiplying a square integer matrix A with det A = +1 (a.k.a. unimodular matrix) [35,36], 
ie. any two frequency vectors w1 and wz can be transformed to each other by such a matrix wọ = Aw ,. Note that the 
transform of A does not change the commensurability of w. 


202211.00236v3 


chinaXiv 


INVARIANT MANIFOLD GROWTH 


Next is the more complicated dummy X-cycle model as shown in Fig. 2, which is acquired by 
twisting the cycle of the first model, i.e. let the R, Z coordinates of the cycle rely on $, instead 
of being constants. In these two analytic examples, we exhibit the technique to construct a field 
by the expected trajectories. At last, a time slice of the magnetic field on EAST is taken as a 
real-world example as shown in Fig. 3 and Fig. 5, with RMP as the non-axisymmetric factor. 


4.1 Non-Möbiusian/Möbiusian saddle cycles 


In this section, a saddle cycle model is constructed as an analytic example. Let By :— RoByo/R to 
simulate a tokamak, where B, denotes the toroidal field magnitude at the axis Ro. Suppose By 
is positive. The trajectories on the invariant manifolds of a Móbiusian cycle is expected to be like 


LA, pines n CH d Ro 
sin 6,,(¢) Zo 

D |e/2mr cos HA) Ro 
i sinds(p)| |Zo 


, (for unstable trajectories) 


X poi (%) E (19) 


, (for stable trajectories) 


where Au, As (independent of p) denote the two eigenvalues of DP", 0,(4), 8s(p) denote the 
corresponding two eigenvectors. If both 6,, and 0, satisfy 0($ + 27) = 0($) + (2k +1)m, k € Z, 
then the cycle is Móbiusian; if they satisfy 0(6 + 277) = 0($) + 2kn, k € Z, then the cycle is 
non-Möbiusian. 


Figure 1: (a, b) show the same Möbiusian saddle cycle, of which the field is constructed by Eq. (26) with parameters: 
(Ro, Zo) = (1.0,0.0), Bgo = 2.5, Ou(p) = p/2 + 7/2, 6(9) = P/2, Aujs = —e*1/5, The sprouts of two 
invariant branches are plotted (red for W", blue for W°). (a) and (b) draw the manifolds on p € [0,27] 
and [0,47] for a half and full poloidal turn, respectively. A trajectory on the unstable branch is drawn for 
three toroidal turns, i.e. @ € [0,67]. 
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Let AX poi(P) := Xpo1(%) — [Ro, Zo]7. For the unstable trajectories, AX oi) is 


S AX (f) = Su Eer p | sells 2m ( eosa] 


dp sin 6, sin Oy 
IN 91] [cosg Inu] gd. 
-]A, [P7 Ka oe asi E zii vil AX pol($) 
0, 2mn u u 2mn 
In |A| o 
= | e m | 9X [y ECH (20) 
2mn u 


It is similar for stable trajectories, with Au, 0,($) replaced by As, Bei), Recall the FLT Eq. (1) in 
cylindrical coordinates is AX ^u) = RByj / Bg (X po1(%), $). Expand RB yo) /By around the cycle, 


R- Ro 


ARB yo /B 
[RBR/ Be ee | ) pot / 2 (6) E H a + O(|AX sall" 


R,Z,$) = | TSAAR 
Seel 9) RcBzo/ Bqo a(R, Z) 
— 
vanishes in this case 


d 
The FLT equation turns out to be AX oi (P) E a AX poi (h) + O(|AX por’). To construct a 


field such that the nearby trajectories are repelled (resp. attracted) at 0, (resp. 05) as expected in 
Eq. (19), it is required that 


In [Au] —6' 
ARB 4 / Bg 2 TW AX poi, if AXy in direction 6, 
"ap zy Xp + O(|AX pil?) = 4p "7 (21) 
a(R, Z) a 7 
2 In| M AXpo1, if AX po) in direction 6; 
0 2mn 


Eigen decomposition is employed to make it possible to let the first orders of the eguations above 
equal at 6,, and 05, respectively. As the cost of such decomposition, 6/, and 6; have to be equal, 
denoted by 6’ afterwards. 


ARByoi / Bo 


s E =e 
ann a E | (22) 


=j 
cos @, cos l | —6' | 
Ee 0, sin, T E Q3) 


In|A 
8 E Oy cos J Du 
sin ĝ, sind, 


Ta fasl 
2mm 


For the left-hand side of Eq. (22), 


Q4) 


IRByoi/Bp — 1 Ea — RBRORBy 9z(RBr)By — Seed 
8.2 ^ B2 ldR(RBz)By — RBzORBy 9z(RBz)By — RBzəzBo ` 


B, has been defined to be RoBgo / R to simulate a tokamak, so ORBy = —Byp / R. Therefore, the 
equation above is simplified to 


Q5) 


IRB poi / Bo u 1 ES + RogBg E 2 E: A R e | 
Q(R,Z) Bp 2Bz--ROgBz RdzBz] By [Bz 0| By lõrBz OzBzl' 


The trajectory Eq. (19) have been preset by known Au, As, 0,($), @s(~), which means the right- 
hand side of Eq. (22) is what we can control. Based on that, the left-hand side of Eq. (22) can be 
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calculated, i.e. the field we want to construct, after which the linearized Bg, Bz field is acquired 
(note Br, Bz at the axis have been preset to zero). 


Se 8 p 2] Pa 2 
ps —lo&Bz zB Z= Za] ^ ONAXroil”) 


E: —0']\ [R- R 
60 0 2 
EC KAN T le J E - A + OAR pol"). ES 
The construction is finished. An example of Móbiusian saddle cycle is shown in Fig. 1. The non- 
Móbiusian case is not shown, because it is easy enough to up Some readers may wonder 


ak RB poi /B $ 


ak 
the higher order terms k > 2, which are determined by IR, ZF = 0 in the left-hand 


Bpo 
a(R, Suo 
side of Eq. (21). 
42 Dummy X-cycle 


In this subsection, Ro, Zo in the last model are replaced (see the trajectory Eq. (19)) with 
R-(¢), Zc($), two functions dependent on 6, respectively. Here, an example expression of R-(¢), Zc(P) 
is given, 


R-(~) =Reii cos (14 a 00) + Rax, (27) 

Zc(%) =Ze sin (up + 09) + Zax, (28) 

where the subscript c is short for cycle, ell for elliptic, and ax for axis. Similar to the first model, 
set Bg(R, Z, $) = Bg(R) = RaxBp.ax/R to simulate a tokamak. 

The zeroth order of field expansion are denoted by Bre, Bzc, Boc, to highlight that they are 


evaluated on the cycle, eg. Bgc = Bg(Re(f), Zc($), $). According to the spiral cycle equation 
defined above, the following equalities hold, 


Xr =R-Bre/Boc =—!Ret sin (up + 00) = RE) (29) 
Xz =R.Bzc/Byc = IZg COS (up + 00) = Z(6). (30) 


The first order of expansion, 0B po) /9(R, Z), is calculated by Eq. (25). 


ap zi d E Be (P PBR Be el) (31) 


AR, Z) Läb: dzBz a(R, Z) 2Bz/Bp 0 
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Expand the poloidal field around the cycle up to first order, 


| 


in Fig. 2. 


Bg 
Bz 


| 


| (6) + [BE] (Zig) + (AX poll) 
[B C OB ol R — R-() 
BE] + se Roy, 24) [2 — 209) + OJA pa?) 


ex] exl 77] - Paterno aD E] + ais, 
ED 


2R//R, 0 
27) /Re 0 


(32) 


The construction is finished. A dummy non-Móbiusian X-cycle with q = m/n = 3/1 is shown 


(a) 


Wy) sprout 
Wy) sprout <—— 


Figure 2: (a, b, c) show the same dummy X-cycle, of which the field is constructed by Eq. (32) with parameters: 


(Rax, Zax) = (1.0,0.0), Boax = 2.5, 1 = n/m = 1/3, (Ren, Zen) = (0.3,0.5), 6,/,($) = i$ +0 = 
¢/34+7/24 7/9, Aus = e*1/5, (a) shows it from the top view, (b, c) show it from the other view. (a, 
b) draw arrows for the two eigenvectors of DP*(4) and their opposite (blue for A < 1, red for A > 1), 
which is acquired by DP~" evolution formula and shall be consistent with 0, Js as designed. The sprouts 
of four invariant branches are plotted (orange for V", blue for W°). In (a, b), the manifolds are not shown 
on $ € [6n—2/3n,6n), in order not to shelter the eigenvector arrows. A transparent torus with the 
corresponding elliptic section is drawn for reference. 
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43 Areal-world example 


The equilibrium field of EAST #103950 shot at 3500ms given by EFIT is taken as background, 
superimposed with a non-axisymmetric field induced by the RMP coils running in n = 1 mode. 
The plasma response is not considered here for simplicity. B, of this shot is negative everywhere, 
and Bpo; at R-Z cross-sections is clockwise. 

On locating the periodic points of the Poincaré map, the simplest discrete Newton method x;,1 = 


xj—h [DG (x;)] ge G (x;) is employed [37], where G(x) := F(x) — I and I is the identity map. To 
locate m-periodic orbits of the Poincaré map P (4) at the section $, our map F is chosen to be the 
mth iterate of P(¢), that is P" (). 

After locating the cycle, one needs to know the DP” (0) of this cycle at every section, which is 
the job of DP*" evolution formula. This formula is a traditional matrix ODE system. The Python 
function scipy.integrate.solve ivp and the Julia package DifferentialEquations.jl have 
prepared a lot of numeric algorithms for such ODEs. Next is to eigen decompose DP" (p), which 
can be done by, for example, the Python function scipy.linalg.eig. The two eigenvectors of 
D'P"'($) and their opposite are the directions towards which the invariant manifolds grow at the 
beginning. 

Recall the invariant manifold growth formula, 


aen: ` E ox) Ié 
às By, Aë J/ 7 


only requires two variables on the right-hand side, RB; /By and ax"/s /ap. In the numeric 
implementation, RB pot / Bọ is linearly interpolated on a regular grid of shape [n pn zn]. For 


Mo, (16 revisited) 


9X''/5 /ðp, different people have different ways to handle it numerically. In the following two 
subsections, two approaches to grow manifolds are exhibited. 


(a) 


Figure 3: (a) 3D visualization of the invariant manifolds of the lower X-cycle of EAST #103950 shot at 3500 ms (EFIT 
+ vacuum RMP). The first wall is drawn as a transparent grey surface, with the top removed in order not 
to shelter the manifolds. (b) Enlarged view of (a) at p = 0 near the lower X-point. 
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43.1 Naive field line tracing 


The simpler scheme is to distribute a line of Poincaré seed points along the eigenvector of 
a periodic point. To compute the arc length s of a branch of W"/*(P"', xo) requires the FLT 
Poincaré orbits used to construct the manifold be ordered, which is achieved with the assistance 
of the DP" ($) eigenvalues in this paper. 


" Expected A 
i W Around a small neighborhood of xo 
X0 
o X 4-99999— — — 
P O XN Oon 
oco xjin the middle 
XQ Not expected yy" P (x1) — xo & Au: (1 — xo) 
To keep the sequence ordered, i.e. |P(x1) — x9| > |xn — xol, 
P one simply needs to let |xy — xo| < Au|x1 — xol. 


Figure 4: Illustration of the numeric algorithm to grow invariant manifolds. 


Suppose xo is a saddle fixed point of the 2D Poincaré map P at an R-Z section. Denote the 
unstable eigenvalue and eigenvector of DP (xo) by Au and ou, If xp is instead an m-periodic saddle 
point, one simply needs to substitute DP” (xo) for D'P(xo). Seed a line of points (x1,..., xy) along 
Uu, equally spaced, i.e. x; 44 — x; is a constant for 0 <i < N — 1. If Ay € N, P(x1) probably falls 


behind xy, ie. the s of P(x1) is smaller than that of xy. This makes it difficult to compute s, 
because the order of s is not certain. Let X be a sequence defined by 


X = (xi,.. XN) SC (P(x1),. ..,P(xNn)) FUN etg 
= (m, ... xv, P(x), ... ,P(xn), P2(x1), ... , P^ (xw), ...). (33) 


It is expected that the s of (x1) be greater than that of xy, the s of P?(x1) be greater than that 
of P(xw), and so on, i.e. the s of P**1(x4) be greater than that of P* (xw). In fact, as long as the 
s of (x1) is greater than that of xy, the conditions for k > 1 are naturally satisfied. Now that 
P(x) ~ xo + Au(x — xo) in the neighborhood of xp if x is in the direction of v,, one simply needs 
to put xy closer to xo than (x1), so that s(P*(x7)) < s(P*1(x4)). By virtue of this fact, the orbits 
are untangled with ease and one gets rid of tentative manifold growth methods [11-15], which 
need to decrease growth step when the local manifold curvature is large. In other words, it is 
ensured that the sequence 


s(X) = (s(x), ... ,8(xn), s(P(x1)), ... ,s(P(xw)), ...) (34) 


is a strictly increasing sequence, which guarantees that it is safe to compute s by simply accu- 
mulating the lengths of segments. Apart from this advantage, this technique to grow manifolds 
applies to hyperbolic fixed points of n-dim maps, n € IN*. 
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The invariant manifolds in Fig. 3 and Fig. 5 are grown by this naive FLT technique, of 
which the computed arc lengths s are expressed by the varying color to let it be more easy- 
to-understand. One can immediately observe the confinement of this equilibrium relies mostly 
on the invariant manifolds of the lower X-cycle jo, although it also has one X-cycle at top. It is 
also known as the disconnected double-null configuration. 

The blue arrows in Fig. 5(c) are drawn by the invariant manifold growth formula, which takes 
0X'/5 [9 and gives 0X'/5 /ds. Evidently, OX? /ds is the growth direction of a manifold. The first 
order central scheme is employed to calculate 9X"/5 /ds (s, 4) by X"/*(s, p) at the two neighboring 
R-Z sections at p +e, 


ax"/s | X'I*(s, o +e) — X"^(s, — e) 


ap 09^ x (35) 


4.3.2 Discretize the invariant manifold growth PDE to an ODE system 


The other numeric way is to transform the invariant manifold growth formula, a partial differential 
equation (PDE) including 0; and Op, to a system of ODEs with s as the evolution parameter. 
Transect the invariant manifold by N R-Z cross-sections at (CARA ; to discretize it, where (CARA , is 
an arithmetic sequence with a common difference Af = $;,4 — $; ranging from 0 to 2m7 — Ad. 
The R and Z coordinates of the manifold at each section contribute two variables of the system, 


x "ie, $;) and xul "Le, $;), totally 2N variables contributed by N sections. Thereafter, Ped *(s, fi) 


and x "(s,$;) become univariate functions dependent on s. The manifold growth formula is 
thereby reduced to a 2N-dim ODE system. 


dX"/s(s, $1) RB yj AX" Í*(s, 61) + IL Ss | 
ds Bo Bn N " 
ep u/s 
dX aS $2)... (Se = JI i + [f (16 discretized) 
dX"/*(s, py) RB yo AX*/5(s, py) + IL e | 
ds By Ap i = 


where AXY/S(s, p;)/ Ap denotes a numeric alternative of JXY/*(s, p)/9p. For example, it can be 
defined as the first- or second-order central difference schemes, 


AX"/S(5, Al. X"/S(S, iss) - X" (s, pii) 


Ap 2Ap í ES 
Ais, éi. X"P(s,di1) - 2X" (s, pi) + XS, pi) (37) 
Ap. Ag? | 


This discretizing method works fine for the invariant manifolds of the outmost saddle cycle(s), 
but not so for the g = m/n = 3/1 island chain, in which case two stable and unstable branches 
grown are almost identical and can not be distinguished. The primary reason is suspected to be 
that the two eigenvectors are so close to each other that the grid granularity is not fine enough 
to offer the difference needed to calculate a correct branch. Although this scheme suffers such 
a numeric instability, the authors still think it is worth introducing, because this scheme directly 
utilizes the invariant manifold growth formula and can be useful for those cases that the two eigen- 
vectors of DP” are not so close. 
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43.3 Denotations and notions explained with the aid of figures 


A periodic orbit for a map is best seen as a whole to reflect the intrinsic homoclinic/heteroclinic 
structure. Then, one can use W” (P, (x1, x2, x3)) to denote all the unstable branches of the X- 
cycle of the q = m/n = 3/1 island chain in Fig. 5. This X-cycle has three strike points (x1, x2, x3 } 
through $ — 0 section, each of which has two stable and two unstable branches. Totally there 
are 6 branches of W° (P, (x1, x2, x3]) and 6 branches of W" (P, (x4, x», x34). One can also spec- 
ify a branch in a more fine-grained way by replacing P with P", (x1,x2,x3) with xj, that is 
w«/s (P3, xj), which represents the two unstable/stable branches belonging to the specific point 
xj, respectively. Obviously, 


W"!5 (P, {xu 22,3) = U W"5(P3x). (38) 
ic(12,3) 


The 2D manifold W"/*(B, y) consists of all the corresponding 1D manifolds yw SP" ($), x(e)) 
for the Poincaré maps P” (9) at all sections d, i.e. 


W"^B,)- LJ (AR,Z,9)E M : (R,Z) E WS(P"(4), x(9))} , (39) 
pe[0,2mz7t) 


where x(p) denotes the (R, Z) coordinates of the cycle y at $ angle, and P($) : R* x R > 
R+ x R denotes the Poincaré map at ¢. where x($) is a 67t-periodic function representing the 
(R, Z) coordinates of this X-cycle y. P(x0) = xi, P(x1) = x», P(x2) = xs. Then, x(0) = x1, 
x(—27)) = x2, x(—47t) = x3. (By is negative in this shot.) Viewing the 3D Fig. 3 and 2D Fig. 5 to- 
gether can help understand the relationship between 2D manifolds W!/S(B, y) and 1D manifolds 
W"Ís(P" (9), x(9)). 

In Fig. 5 (b), some regions enclosed by the invariant manifolds of Jj are filled by colored 
scatter points. Points are marked with distinctive colors, and then mapped by P~!. How these 
regions are connected by field lines is reflected through the color pattern of scatter points. Let x79; 
be the strike point of "mee crossing the $ = 0 section. All the points of transversal intersection 
of the two manifolds W"(P,xj,,,) and W*(P, xiow) (the red and blue curves, respectively) are 
homoclinic to Sins, 

Although the fluxes in the colored regions filled by scatter points are the same, e.g. the @ and 
@ regions in Fig. 5 (b), this flux value probably differs from that of the uncolored regions, e.g. (3) 
and (4). One should always be careful of this fact. 
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^ 
d e 
d e 


-0.4 


-0.5 


1:6 1.8 2.0 2.2 24 


WP, x) 


W(P, x) 7 


W(P, x) 


Xı 
x 


a? Px) 


1.623 


1.626 1.627 


Dense seed points of FLT are 
put here, and then mapped by 
P ~! once and once again. 


1.8 1.9 2.0 


Figure 5: (a) Poincaré plot at p = 0 of EAST #103950 shot at 3500 ms (EFIT + vacuum RMP). Some invariant 


manifolds are grown and plotted. (b) Enlarged view of (a) near the lower X-point. Dense scatter points with 
color are presented to show how the regions are connected by field lines. (c) Enlarged view of (b) near xy. 


The blue arrows are drawn according to Eg. (35). 
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5. DISCUSSION AND CONCLUSION 


5.1 Comparisons with existing works 


This subsection compares various approaches developed to study invariant manifolds, albeit 
obvious. The comparisons of our D'P" evolution and invariant manifold growth formulas with 
others' works are presented in this subsection subsequently. 

Most of existing research has been satisfied with the Floquet's normal form. Floquet theory is so 
successful that people cease further exploration. Floquet theorem governs how DX pol (Ps, Pe) varies 
with Ge and guides how to solve the system in a simpler way. However, it remains unclear how 
DX pol(Ps, Ps + 2mn) changes with ps, which is essential to lessen the computational resources 
needed to solve for the initial growth directions of stable and unstable manifolds. The most 
similar work to DP" evolution formula is given by Tsutumi [38]. Tsutumi considered the following 


linear system: 

= = A(p, t)x, —oo < p,t < +00, (1.1) in [38] 
where x = x(p,t) is a complex n-column vector and A(G,t) is a complex n x n matrix function 
T-periodic in $. Theorem 1. in [58] is repeated as below: 


There exists a monodromy matrix of (1.1) which does not depend on t ^ if and only if there exists a 
matrix function T(G, t) which is defined on co < p,t < œ, T-periodic in $, and satisfies 


d d 
FADO- a, (0 + [A(6, D, FG, 0] = 0 (1.2) in [38] 


Tsutsumi's Eq. (1.2) above is more complicated than our DP” evolution formula (11) because 
A(¢, t) relies on both $ and t. If JA /0t vanishes, one can easily observe that T is governed by the 
same equation as DP*" evolution formula. Tsutsumi did not explain in [38] what T can be and 
how to construct it. His attention was paid to the condition of A($, t) such that the monodromy 
matrix does not depend on t. This condition is revealed in his theorem by the existence of the 
matrix function I which satisfies Eq. (1.2). 

The works similar to our invariant manifold growth formula are collated into Table 2, in which 
the work of S. S. Abdullaev is discussed in the next paragraph and Appendix B in more detail. 
As shown in Table 2, the classical invariance equation (1.7) in [10] for n-dim maps might be too 
general for fusion scientists to be useful. T. E. Evans, J. P. England, and J. M. Ottino et al. focused 
on numerical or experimental methods rather than analysis. 

Abdullaev [30] has made some contributions to establishing the formula of stable and unstable 
manifolds. Here we point out the following facts worth noting about his work: (a) His theory 
relies on the assumption that the perturbation to an axisymmetric magnetic field is small so that 
the FLT system can be transformed to a perturbed Hamiltonian form; (b) The Poincaré integral 
is employed to acquire a first order approximation to Poincaré map. By contrast, our Poincaré 
map is exact; (c) The final formula F Gig, x,y) = 0 is implicit instead of explicit; (d) The (x, y) 
coordinates appearing in the formula are not the standard (R, Z) cylindrical coordinates. The x- 
and y-axes are even probably not perpendicular. A more detailed comparison equipped with 
equations is put in Appendix B. 


“In [38], "there exists a monodromy matrix of (1.1) which does not depend on t” is equivalent to that “the internal struc- 
ture of every monodromy matrix of (1.1) does not depend on t”. By internal structure, Tsutumi means the characteristic 
multipliers of a matrix and their algebraic and geometric multiplicities. 
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Manifold expression comparison 


Table 2 
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5.2 Conclusion 


| Beginning 


Attempt to locate X/O- | 


| points at a R-Z section 


By classical Newton method or 
any other methods you favor 


What is the geometric 
meaning of |DP+™| ? 


` X/O-points located at the R-Z section, 


+m 
Calculate DP Nu and DP” ( D)lp=o acquired | 
The eigenvalue of 
How? DP+t™() determines 
The invariant manifolds grow from WEE manifold 
: e F m 
Ee eun LED nuc (9. z| is unstable/center/stable. 
No S 7 | 
Simply integrate 
ODX po1/0 wrt. A Yes 
for 2mm 
d RBy 0 xu/s 
io 77) = [A(p), DP $^ ($)] ax"s Bo op 
ðs £l RByoi 8 axu/s I 
Ka Bo 0b 2 


Figure 6: Thought process map of this paper, where the contents in grey boxes are not new. 


Although the magnetic fields in tokamaks are mostly considered axisymmetric, almost all 
kinds of auxiliary heating schemes, except the ohmic heating of the central solenoid, are strongly 
localized and non-axisymmetric. The topology of magnetic field plays a dominant role in the 
behaviour of the confined plasma and the scrape-off layer. Motivated by the curiosity about the 
intrinsic characteristics of general 3D vector fields (not necessarily divergence-free), an analytic 
theory on the invariant manifolds of cycles is established in this paper, where the short-period 
cycles are regarded the skeleton of fields. 

The invariant manifolds of saddle and parabolic cycles grow from the eigenvectors of DP ^" (q). 
By D'P-" evolution formula (11), one gets rid of repetive 2m7t-integration of Eq. (3) for D'P^"' (9) 
at every angle $. The primitive FLT ODE system (1) is extended to the invariant manifold growth 
formula in cylindrical coordinates (16) by analyzing the relevant differentials. 

It is proposed to use the notion of invariant manifolds of outmost saddle cycle(s) to characterize 
the magnetic topology at the plasma edge when the field is strongly non-axisymmetric. If a 
tokamak operates in the single-null (resp. double-null, or somehow more strange) mode, there 
is one (resp. two, or more) outmost saddle cycle(s). The transversely intersecting manifolds are 
suspected of causing the spiral ribbon-like pattern of heat deposition on the divertor found in 
tokamak experiments, which needs further verification. In the scrape-off layer, how to disperse 
the particle and heat fluxes before they land on the divertor is also an interesting problem, worthy 
of more investigation. With the drift effects taken into account in the future, the heat load pattern 
observed by diagnostics like infrared thermography is supposed to be more consistent with the 
divertor regions covered by the invariant manifolds of outmost saddle cycles. For both tokamak 
and stellarator communities, the transport issue at the plasma boundary is essential to control 
the heat load below the material limit. 
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APPENDIX A. PROOFS 


A.1 The geometric meaning of [Dx pol(%s, Pe) 


ge R(V -B Bolo, D 
[PX us, fe)| = exp (I ( ) ao ele. (8 revisited) 
a Bo WI 
Proof. 
$e Pe 
[PX poi elg Ad opi / On(RBg/ Bg) + 94 (RBz [By)äp) (40) 
Pe 1 
—exp d zz (RBpOnBr + 92B2) + BrBp — R(Brdr + Bz0z)B4) dé (41) 
* 74 
since V - B = ðrBr + 0zBz + 1/R (Bg + 0$ Bal (42) 
be 1 
= eel JA go (RBp(V - B) — BydpBy — R(Brðr + Bz9z)Bp) Z (43) 
Jos Be 
To further simplify the expression above, prepare the following differentials. 
A X , 1 
VBo =bRORBy km e70z Bg EE ep R99B9 (44) 
dl —égdR + ezdZ + égRdp (45) 
VByp - dl =dgBgdR + de Bad? + da Bod (46) 
VB,:dl/d$ =dpBy Xr +0zBp Xz +əB 47 
p: dl/dp =ðrBp Xr zBy Xz 9 Bp (47) 
=RBr/By =RBz / By 
Reorganize the terms at the right-hand side of Eq. (43) in a much more compact way, 
be R(V.B) 1 RBR RBz 
DX pi zi E (^ B + Rap oe T an. ) do (48) 
| po Jos By By gap By $ By d 
be R(V-B) VBg-dl 
= exp / — dọ (49) 
| «Be Bodo 
% R(V.B te VB, - dl te R(V.B de 
=ep| f VD ap NEG: -eo( / s aen ) (50) 
% Bo $ Be Jos By $s 
—~S 
=V In By-dl 
(Here In is the complex logarithm instead of the real one. So it can handle negative numbers.) 
Pe R(V-B Bolos 
= exp (/, x ay) Bele. (51) 
Jo Bg Bola 
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A.2 DX to DP formula 


1 RBR 
=a : 
C7! DX C; = DX pa | . (10 revisited) 
1 —RE d 
Bo J | eng 
cos $ —Rsing 
Cs and C, are the matrix | sing R cos p | evaluated at the starting and ending points, respec- 
1 
tively. 
Proof. Consider the differential relationship in cylindrical coordinates, 
dx =  (dR)cos$ — R sin dọ dx cos $ —Rsin$| [dR 
dy= (dR)sinp +Rcosfdp = dy| = |sing Reos@ | |dZ (52) 
dz = dZ dz 1 dé 


Suppose the solution X(xo, t) crosses the R-Z cross-section e at the time T. The relevant differentials are related by 
the equation dX (xo, T) = DX(xo, T) : dxo as detailed below 


dX —DX dxo, 
dX; | dxox 
dX,| =DX |dxo, ^ , 
dX; | dxoz 
cos p —Rsing dXg] cos p —Rsing dxor 
sing Rcosp dXz| =DX |sing Rcos$ dxoz|, 
1 end dX 4 1 start dxog 
dXg] dxor 
C, |dXz| =DX C, |dxoz| , (53) 
dX, | dxop 


where xo and X(xo, T) denote the starting and ending point, respectively. A similar differential analysis is done for 
X poi Qs, $e) as below, 


[ RB 
mei (54) 
pol dxoz dXz = Bax, 
? end 
1 — A dër 
E i = RB; dXz (55) 
L By end dX 
; 1 — RBR dxor 
FOR M Ei | C;' DX C, |dxoz |, (56) 
L Bo end dxop 
The equation above always holds for any dxor, dxoz, and dro, hence the corresponding coefficients shall equal, 
1 — RBR 
* Bp 1 
[px pol | = C; DX Cs. (10 revisited) 
P 1 — RBz 
By end 
A.3 DP+" evolution formula 
EE TJ " 
a (4) = [A($), DPH" (O)|, (11 revisited) 


Proof. 
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Integrating Eq. (3) w.7.t. pe from $ to ps; + Ad, 


I rpst+ag a 
DX, ds + AP) — DX ae) = L ag DX) dhe (57) 


m 


—M— 
=A(fe) DX poi (Pse) by Eq. (3) 


Differentiating both sides of Eq. (57) w.r.t. os, 


ag, DX rab $s + A) =A (Ps + AP) DX por (Ps, Ps + Ap) — Als) (58) 
ds +Ap A 
dE HEH e 


d d 
Though it is known how to calculate ae P pol (Ps, Pe) as shown in Eq. (3), little is known about a pol (Ps, Pe). 


Notice that DX pol (Ps, Pe) is the inverse matrix of DX pol (Pe, gs). The following Eq. (62) from [39] is useful to solve for 
(KY by K(x) and its derivative K'(x), where K(x) is a univariate matrix function. 


(KK) =' =0 (60) 
=K'K=' + K (K=!) (61) 
> (KT) = kk (62) 


9 
Borrow the Eq. (62) to solve for n a pol (Ps, Pe), as shown below 


d d E Eg. (62) d os 
ag, PX ros Ai = 59, PEro Pep) TU (DX Per al 35; DX ot fe #3) (DX (her s)) 
—— 
A (s)DX poi (ers) 
= — (DX (60,95) Alps) DX poi (Per PDX (9e, Ps)) 
xx DX (ds, Pe) A(¢s). (63) 


Next, calculate the total derivative of DX po (ps, Ps + Aq) wr.t. ps. Let pe = Ps + AQ. 


CS A(de) DX pol (ps Pe) A(gs) 


eo 
DX pil srs + AG) CAQs + A9) DX pilots + AD — AG) + [^ So [A9 DX (fs 9] dde 
dé, "UU i Jø dé m 


—A($s + AG) DX poi (Ps, ds + A) — Als) — (DX poi (hs, Ps + Ad) — I) A(s) 
—A($s T Ad) DX poi (Qs, Ds + Ad) m DX poi (Ps, Ps + AP)A($5) (64) 


For trajectories of m toroidal turn(s), A(s) and A(Gs + 2m7) are identical since A is a periodic function on the cycle. 


A(4s) = A(ps + 2mm) (65) 


Let Ap be 32m, then DX il, Ps + ^d) can be substituted for DP" (ps). Eq. (64) is simplified to 


d 
d$; 


DP" (ps) = Alps) DP^" — DP=™ A(ps) 


= [A(6s), DP, (11 revisited) 
One might suspect whether Eq. (11) works for DP ^", the inverse of DP”. Consider DP™ = (D'P"')-1. By Eq. (62), 


d —m mi— m d GA keng my— m m my— 
ag P OP) '(pP") (DP)! = —(DP™) ! (a DP" — DP" ADP")! 


= [A($s), DP" (9,)) , (66) 


which indicates that the formula (11) applies for both P” and its inverse P". 


23 


INVARIANT MANIFOLD GROWTH 


AA DP" eigenvector evolution formula 


6' = (|, B VA- DPH" 8 B v) B (DP* 


Em y, 


(13 revisited) 


Proof. Firstly. the eigenvectors are parameterized by 0 as shown below, 


DP*" vj(9) = Af), let o; = IO Alum, ie {1,2} (67) 
1 

which is safe because, according to the A-invariance property, the eigenvectors would not become complex on saddle 

cycles. 


Concatenate the two eigenvectors into a matrix V, then 


M 


DP" Tout, v2(4)| = [A191(4), A292(4)] = V kä D |. 
—— - 


(68) 
V <a 
Differentiating Eq. (68) on $, 


(bpiey V 4L Dp*n y! = VA -+VA' 
V' can be simplified by 


hi vB d 
vV = | E v " (70) 
w 
=0' 
O and A are diagonal and hence commute, therefore Eq. (69) becomes 


(69) 


(DPp+"y' V —V'A — DPH". N! VA! 


= F =] VAO’ — DPH 


(71) 
em l m VO + VA. 


(72) 
For cycles, (DP”)’ = LA, D'P"] implies that A = 0, A5 = 0 as shown by the Eq. (15), then Eq. (72) is simplified to be 
(Dpimy X 


lk paeem. poe 


(73) 
A.5 Invariant manifold growth formula 
RB poi ox"/s 
u/s B, 39 
SS $ ? (16 revisited) 
de RBpol axu/s 
"| By oo 


2 
Proof. The relationship among the differentials ds, d$, dR, and dZ is shown in Fig. 7. To solve for ds/d(, the 
differential ds is given as below, 


2 u/s RBR u/s ` u/s RBz u/s ` 
ds” = | Xs Guru AR Ae ($--d$)] + (X2 Cpt -py E. (s,ġ +dọ) (74) 
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X"/5(g -- ds, $ + d$) 


Figure 7: Geometric diagram to show the relationship among the differentials ds, Ap, AR, and dZ. It is supposed that 
there exists a hyperbolic cycle at bottom, from which an invariant manifold grows. (a) shows the differential 


involved in FLT, (b) shows an iso-s curve and an iso-$ curve of the manifold. 


Dividing both sides by d4?, 


( 


dé dé Bo dé 


2 2 
Ep, AN" (RE gu 
B — 36 By — 3$ 
2 2 
ds | |(RB& 9X (pg 9x7 
a ! | B, m By ap A (the sguare removed) 


Express the differential x *(s +. ds, $ + d$) — xev *(s, $) with ds and d$, 


=RBr/Bpdy (simply along the field line, see Eqs. (1) ) 

EE /s ax! /s 
R ds + —R do 

de oo 


(divide both sides by df) 


Xil + dsp + dp) — XW/(5,9) = 


RBR ` ON ds E 
Bp — s dp ^ dp 
WY 


Eq. (77) 


By ð$ aah 
u/s 2 u/s 2 ds ` 
RBg _ 9Xg , ( RBz  9Xz 
e Bo Er] N By Er] 


The Z component as * /ds is similar. 


u/s u/s 2 u/s u/s 
dee E (5,9) -XR (6,9 + dd) | EM) " E (6,9) — X2 (6,9 +44) | RBz 


By 


(75) 


(76) 


(77) 


(78) 


(79) 


(80) 
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APPENDIX B. COMPARISON WITH ÄBDULLAEV'S WORK 


Abdullaev firstly transforms the standard cylindrical coordinates (R, Z) to the canonically conju- 


gated variables (z, pz): 
Z As RAg l 
= = 1 
Ry’ VBR Ty BoR2 Caley 


where By is the toroidal magnetic field strength at the major radius Ro, and A = (Ag = 0, Az, Ag) 
denotes the vector potential of magnetic field. The normalized poloidal flux pp serves as the 
Hamiltonian. However, this transform is not always viable when the magnetic field is non- 
axisymmetric (the definitoin of pp requires a closed flux surface diffeomorphic to T?), which 
limits the application of Abdullaev's method to those near-integrable systems. After that, the 
FLT equations are presented in the following Hamiltonian form: 


dz m Opp dpz u dpp 
dp ðpz dé dz 


(2) in [30] 


A non-axisymmetric magnetic perturbation changes the original axisymmetic poloidal flux ps 
to: 

bp = VO (z, pz) tel (z, pz), (6) in [30] 
where € is a dimensionless perturbation parameter dictating the relative strength of perturbation, 
and py (z, pz, $) denotes the perturbation magnetic flux. Such a denotation implies an additional 
numeric step to transform the perturbation field Bper:(R, Z, $) to ps? (z, pz, $). By comparison, 
B = Bo(R, Z) + Bpert(R,Z,4) is taken as a whole in our theory, without need to distinguish the 
perturbation field from the field to be perturbed. 

The Poincaré map in [30] is defined to be a full poloidal turn map, with the Poincaré section Xs 
consisting of two segments of the C- and y-axes. (č, 4) is a rectangualr coordinate system centered 
at the X-point (see Fig. 1 in [30] ). Note that y, varies along the ¢-axis. Under the perturbion ps, 
the Poincaré map (Pr, Yk) — ($i. i1, Yk+1) has the following classical general form: 


.QOP (dy + 79 (x) ; kai) 


Pk+1 =Y] m 
aP (p£ (He) Yas) irs 
Pk+1 =k + € = E TEEL E me [a (y) +4 të 
k+1 


where the upper and lower signs correspond to the $-increasing and $-decreasing Poincaré maps, 
respectively. P($; wv) is an integral of e) along the unperturbed trajectory, a.k.a. the Poincaré 
integral, 


Rou 
P) | 99 (2 (6:8). P: (9:4),4 +4) dd, (9) in [30] 
—74(%) 


Be careful that the equality of Eqs. (8) only holds when e is infinitesimal, otherwise it is suggested 
to use z instead of — in Eqs. (8) or add a remainder to omit the high order terms for correctness. 
Eqs. (8) are the basis for many other equations in [30], which probably also needs to change to ~ 
when applied to a realistic e. 

Let Yk+1 = 0 at $y,1 — +o in Eqs. (8), then Abdullaev immediately obtains the following 
implicit form of the subset of these two manifolds on Xs: 


FS, y) = YF ea, P + zq();0) = 0, (59, 60) in [30] 
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where the upper and lower signs correspond to the stable and unstable manifolds, respectively. 

For the points of these two manifolds not on X, they are mapped to 2, along the field lines 
by Eqs. (36) in [30]. Under the assumption that the magnetic perturbation at the X-point is small 
and can be neglected, ay, ay coefficient in Eqs. (36) are omitted to obtain the following simplified 
implicit manifold expression: 


: d 1 Q 1 x A 
FO) (p, x,y) = yxy F STM + —In In ;0) =0, (70) in [30] 
(Pay) = vr teg FPE Se dal 2y ly 
Note that (x, y) in [50] are not the standard cylindrical coordinates (R, Z), but defined by 
ag + Br] —a6 + By 
x= j = art of (35) in [30] 
V 2Y y Dë P 


where «, B come from the second order term of the expansion of h in the (€,7,4) coordinate 
system. 


nic, 1, p) = Up (z, Pz, ol — ps (Zs, ps) 
2 2 
& € [ag(g)€ + a(n] — =? + Ca (27) in [30] 


Remind that x- and y-axes are not necessarily perpendicular. The slope of x-axis in (4, Cl plane is 


B/x (y = (—aé + Bi) / /2y/ = 0 on x-axis), while that of y-axis is —B/a (x = (ač + Bi)/ /2y = 0 
on y-axis). The product of their slopes equals — 8?/4?, which does not necessarily equal —1. 


REFERENCES 


[1] D'haeseleer W D et al. 2012 Flux Coordinates and Magnetic Field Structure: A Guide to a Funda- 
mental Tool of Plasma Theory (Berlin: Springer-Verlag) 


[2] Finken K H et al. 2005 The Structure of Magnetic Field in the TEXTOR-DED (Jülich: 
Forschungszentrum Jülich GmbH) 


[3] Liang Y et al. 2007 Phys. Rev. Lett. 98 265004 

[4] Liang Y et al. 2010 Phys. Rev. Lett. 105 065001 

[5] Liang Y et al. 2013 Phys. Rev. Lett. 110 235002 

[6] Cary R C et al. 1983 Ann. Phys. 151 1-34 

[7] Nardon E et al. 2007 J. Nucl. Mater. 363-5 1071-5 

[8] Evans T E 2015 Plasma Phys. Control. Fusion 57 123001 


[9] Wiggins S 1994 Normally hyperbolic invariant manifolds in dynamical systems (New York: 
Springer Science & Business Media) 


[10] Haro A et al. 2016 The Parameterization Method for Invariant Manifolds (Cham: Springer) 
[11] Krauskopf B et al. 1998 Internat. J. Bifur. Chaos 8 3 483-503 
[12] Krauskopf B et al. 1998 J. Comput. Phys. 146 404-19 


27 


INVARIANT MANIFOLD GROWTH 


[13] England J P et al. 2004 SIAM J. Appl. Dyn. Syst. 3 Iss. 2 161-90 
[14] England J P et al. 2005 SIAM J. Appl. Dyn. Syst. 4 Iss. 4 1008-41 
[15] England J P et al. 2007 Internat. J. Bifurcat. Chaos 17 Iss. 3 805-22 
[16] Evans T E et al. 2002 Phys. Plasmas 9 4957 


[17] Kuznetsov Y A 2013 Elements of Applied Bifurcation Theory (New York: Springer Science & 
Business Media) 


[18] Kuznetsov Y A et al. 2019 Numerical Bifurcation Analysis of Maps: From Theory to Software 
(Cambridge: Cambridge University Press) 


[19] Frerichs H et al. 2015 Phys. Plasmas 22 072508 


[20] Frerichs H G et al. 2010 Three-dimensional Plasma Transport in Open Chaotic Magnetic Fields: A 
Computational Assessment for Tokamak Edge layers (Jülich: Forschungszentrum Jülich GmbH) 


[21] Xu S et al. 2018 Nucl. Fusion 58 106008 

[22] Zhou S et al. 2022 Nucl. Fusion 62 106002 

[23] Loarte A et al. 2017 Fusion Eng. Des. 122 256-73 
[24] Frerichs H et al. 2020 Phys. Rev. Lett. 125 155001 
[25] Jia M et al. 2021 Nucl. Fusion 61 106023 

[26] Pitts R A et al. 2019 Nucl. Mater. Energy 20 100696 


[27] Ottino J M et al. 1989 The kinematics of mixing: stretching, chaos, and transport (Cambridge: 
Cambridge university press) 


[28] Roeder R K W et al. 2003 Phys. Plasmas 10 3796-99 
[29] Evans T E et al. 2005 J. Phys.: Conf. Ser. 7 174 
[30] Abdullaev S S 2014 Nucl. Fusion 54 064004 


[31] Jacob Palis Jr et al. Geometric Theory of Dynamical Systems : An Introduction, translated by 
Manning A K (New York : Springer-Verlag) 


[32] Teschl G 2012 Ordinary Differential Equations and Dynamical Systems (American Mathematical 
Society) 


[33] Hirsch M W et al. 2013 Differential Equations, Dynamical Systems, and an Introduction to Chaos 
— 3rd ed. (Elsevier) 


[34] Nelson R B 1976 AIAA J. 14 1201-05 


[35] Broer H W et al. 2010 Chap 6 KAM Theory: Quasi-periodicity in Dynamical Systems. In Handbook 
of Dynamical Systems 249-344 (Elsevier) 


[36] Olikara Z P 2010 Computation of Quasi-periodic Tori in the Circular Restricted Three-body Problem 
(Purdue University) 


28 


202211.00236v3 


chinaXiv 


INVARIANT MANIFOLD GROWTH 


[37] Miller J R et al. 2000 Phys. D: Nonlinear Phenom. 135 195-211 
[38] Tsutsumi M 1975 Proc. Japan Acad. 51 645-64989, 4, 1010-21 


[39] AT 2015 https:/ / math.stackexchange.com/ q/ 1471836 ver 2015-10-09 (online: Mathematics 
Stack Exchange) 


29 


